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Strong approximation for It6 
stochastic differential equations 


M. Namjoo* 


Abstract 


In this paper, a class of semi-implicit two-stage stochastic Runge-Kutta 
methods (SRKs) of strong global order one, with minimum principal er- 
ror constants are given. These methods are applied to solve It6 stochas- 
tic differential equations (SDEs) with a Wiener process. The efficiency of 
this method with respect to explicit two-stage It6 Runge-Kutta methods 
(IRKs), It6 method, Milstien method, semi-implicit and implicit two-stage 
Stratonovich Runge-Kutta methods are demonstrated by presenting some 
numerical results. 


Keywords: Stochastic differential equations; Strong approximation; Runge- 
Kutta methods. 


1 Introduction 


In recent years, a great deal of concern has been raised regarding the study of 
SDEs as an important area of research. Many phenomena in science and en- 
gineering have been modeled by deterministic ordinary differential equations 
(DODEs). However, some of the parameters and initial data are not known 
with complete certainty due to lack of information. Therefore, to represent a 
more accurate model of the behavior of such phenomena they usually should 
be modeled by SDEs. Some areas where SDEs have been used extensively in 
modeling phenomena include chemistry, physics, engineering, mathematical 
biology and finance (see, for example, [5], [7]). Since explicit solutions are 
known only for a few equations, the study of numerical methods have become 
more important and these must be designed to be implemented with a certain 
order of accuracy. Consider the autonomous It6 SDE given by 


dy(t) = go(y(t))dt+ ga(y@))dWt), y(to)=yo, té€ [to,ty],  () 
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where go and gj; are real-valued functions which are called the drift coeffi- 
cient and the diffusion coefficient, respectively, and W(t) is a one-dimensional 
standard Wiener process, whose increment AW (t) = W(t+h) — W(¢) has a 
Gaussian distribution with mean 0 and variance h, i.e. W(t +h) — W(t) ~ 
N(0,h) = Vh N(0,1), and the solution y(t) is an It6 process. A Wiener pro- 
cess (named after N. Wiener) is sometimes called Brownian Motion, which is 
a term used to describe the phenomenon of the erratic behaviour of a particle 
in a liquid, acted on by random impulses, in the absence of friction. Equation 
(1) can also be written as a stochastic integral equation 


y(t) = yo + / go(y(s))ds + / au(y(s))aW(s), 


to to 


where the first integral is a mean square Riemann-Stieltjes integral and the 
second integral is a stochastic integral which can be interpreted in many 
ways (see [10]). The two most studied interpretations are due to It6 and 
Stratonovich that depend on the points of the partitioning in which the inte- 
grand is evaluated. If the lower end point t, is chosen, it leads to It6 integral 
and if midpoint (ty, + tn41)/2 is chosen, it leads to Stratonovich integral. 
The Stratonovich interpretation follows the common rules of integral calcu- 
lus, while the Ito formulation has the advantage of preserving the martingale 
property of Wiener process. It is always possible to switch from one interpre- 
tation to the other, because an It6 SDE can be converted to a Stratonovich 
SDE (and vice versa) by means of the following formula (see [5]) 


Go(y) = go(u) ~ 59 (v)ar(y), 


where equation (1) is in the Stratonovich form when gp is used in place of go. 
There are different numerical methods for solving these kinds of differential 
equations (see, for example, [1], [6], [8]). Numerical methods for SDEs are 
recursive methods where trajectories, in other words, the sample paths of 
solution are computed at discrete time steps. These methods are classified 
to strong and weak. Only strong convergence will be considered in this 
paper. Strong convergence is required, when each trajectory of the numerical 
method must be closed to the exact solution. Formally, if yxy is the numerical 
approximation to y(t) after N steps with constant stepsize h = (t; —to)/N, 
then yy is said to converge strongly to y(ty) with strong order p if there 
exists C' > 0 (independent of h) and 6 > 0 such that 


o( 


An outline of this paper is as follows: In Section 2, the semi-implicit SRKs for 
SDEs are introduced, moreover order conditions for a class of SRKs with order 
one are stated. In particular, the new class of semi-implicit two-stage SRKs 
for SDEs with minimum principal error constants is constructed and the fixed 


yn — y(tn)|) <Ch?, he (0,6). 
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point iteration algorithm will be used to improve the semi-implicit method. 
In Section 3 we have some numerical results which show the efficiency of this 
method. 


2 The semi-implicit Ito Runge-Kutta methods for SDEs 


The most famous numerical method that can be obtained from a stochastic 
Taylor expansion is Milstein method. This method for the SDE problem (1) 
is given by 


1 
Yn+1 = Ym + hgo(Yn) + Sigi(Yn) + 5(r” — h) 91 (Ym) 91 (Yn); 


where Jj = W(t, +h)—W(t,) with h = (t¢—to)/N for some integer N. This 
method converges with strong order one as long as E(y?) < oo, and go, 9, 
91, 9; and gi’ satisfy a uniform Lipschitz condition. Higher order numerical 
methods can be obtained by truncating farther terms of the stochastic Taylor 
expansion. This technique involves considerable complexities in implementa- 
tion because of the approximation of higher order stochastic integrals and 
the evaluation of high order derivatives of both the drift and diffusion coef- 
ficients. Thus, it is important to be able to derive derivatives free numerical 
methods and this leads to SRKs. For the SDE (1) SRKs is given by (see [2]): 


0 1 . 
¥i= tn t+ 22 g0(¥3) + Zi al¥), 1=1,2,..-,8, (2) 
j=l j=l 
Yn+1 = Yn + 2°) go(¥5) F S291 (Yj), 
j=l j=l 


which can be represented in tableau form as 


ZO) | 7Q) 

au | ay 
where Z(*) = ae for i,7 = 1,2,...,s and glk)? — (26, = _, hh) repre- 
sents for k = 0,1. Here Y,,..., Y, represent the internal stage of the method, 
and Yy,41 is the update of the numerical solution at the end of the current 
step. Since (2) is a generalization of the class of Runge-Kutta methods in 
deterministic case, for consistency the stepsize will be included in the param- 
eter matrix associated with the deterministic components, so Z) = hA and 
ZO" = ha’, while Z“ and 2” have elements that are arbitrary random 
variables. In order to derive methods with strong global order one, the ex- 
istence of stochastic Taylor series expansion of the SRK method in the Itd 
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case and the Ito Taylor series expansion of the exact solution is necessary. 
By comparing these two expansions, the local truncation error over one step 
with an exact initial value can be written as (see [3]): 


L(to +h) = y(to +h) —YV(to +h) = Do elt) Ft)yo, 
teT* 


where e(t) and F'(t)yo are called the local truncation error coefficients and the 
elementary differential for tree t, respectively and T* is the set of bi-coloured 
rooted trees . Assuming certain conditions on the cofficients of the method 
and satisfying Lipschitz condition for the drift and diffusion cofficients SDE, a 
method will have strong global convergence of order one if it has strong local 
order one and mean local order one (see [3]). In [9] the order one conditions 
for a class of IRKs in the form 


ZO—hA, 27 =hat, 20 = VaEBO + 7,8, (3) 
T T T 
BO ap age. 


are given, where A, B“) and B®) are s x s real matrices, and a? = 


th ag 


are row s-dimensional vectors. In fact a SRK method of the form (3) will 
have strong global order one if (see [9]) 


ate =1, 
ye =0, 
OM = 1, 


yO? BMe = —1 


2? 


yO)" Be 4+ 727 BVe =0, 


4 
72) Be = 1, @) 
aT BYe = 0, 
vy)" Ae = 0, 


yO" (BMe)? +07 (Be)? + 2y)7 (BM e)(Be) = 0, 


(1)? Ba)? (O* 9Oyr nk (gO poe. FORM) ao 
y ea saad a | = 5 ) 


Here e = (1,...,1) € R® and multiplication of vectors are componentwise. 
If the matrices A, B“) and B®) are strictly lower triangular, then the method 
(3) is said to be explicit, while if A, B™ and B®) are lower triangular, then 
the method (3) is said to be semi-implicit. A family of two-stage explicit 
SRKs of the form (3) with minimum principal error terms can be presented 
by the following tableau (see [9]): 


0 0 0 0 
0 0| -}(VR- A) 0 
[nh O| —-vh Vht Ji 
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which is called ‘E.M1’, and has the principal error constants 


and the other family of two-stage explicit methods satisfying (4) with mini- 
mum principal error constants can be presented by (see [9]): 


0 0 0 0 
0 0| AvAtA) 0 
[rh o| va —Vht+ A 


Also the It6 method (see [2]) that is a derivative free version of the Mil- 
stein method with strong global order one, can be presented by the following 
tableau: 


0 O 0 0 
0 0 Vvh 0 

2 2 
ho Of A-A(AY =D) BCA -1) 


This method is called ‘I R.K’ and has the principal error constants 


In [1] a class of semi-implicit and implicit Stratonovich Runge-Kutta meth- 
ods of strong order one with minimum principal error constants for SDEs is 
constructed. More precisely, this class of semi-implicit and implicit two-stage 
Stratonovich Runge-Kutta methods with minimum principal error cofficients 
can be presented, by the following tableau (see [1]) 


3473 34-73 
ce ee 
3 6 h 3 1 6 Ji 


| gh gh 1 31 
and 
=P h ght si Ji ron 
| sh gh at x1 


which are called ‘S7M1’ and ‘IM’, respectively. In order to generalize the 
above explicit IRKs to semi-implicit case, consider s = 2, hence the matrices 
A, B® and B®) will have the following forms: 
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1 2 
a21 G22)’ os) os) ' oe) oe) 


Now by system equations (4) and according to the structure of matrices A, 


B® and B®) of the above form and by MAPLE, we have the following 
system equations of ten equations with fifteen unknowns: 


Es a 
a?) 4?) = 1, 


1 
Oe) py 4 olD = 1, 


4,0), @),.@) , 2) , Or), (),) , (2) ,0) 

V1 Oy +99 boy +79 b59 + bly + 9° bay +79 O39 = 9,~ 
(2),(2)  (2)2(2) 4 (2),(2) — 1 

aa (ciel 

abo) + a2b6) + a2b%y) =0, 


yO arr + Yaar + 8? a22 59, : (5) 
(),M? . 2) (,@) . p@ (1),(2)? | 1) (2) 4 2) (2) (1) 2 (2) 
nob +75 (01 + oi) +77 bry +75 Cs +oY)) + yy by bry 
2 2 2 1 1 
+2r2? (u + 02) (16 +0) =0, 


(2), (1), (2) (2) (,(1) (2) (1),(2) (1),(2) (2), (1) (2) (1) 
29, Oi bry +79 Cs bry + bgg boy + 2b39'ba9° + b31'by1' + boo b3) ) 

1),(1)? 1 1), 1), 1),(1)? 1), (2)? 
10? 46 (AD + UD) 4 oD? 4 af 

1 2),(2 2), (2 1) (2)? 
a6? (WDD + AD) + 400" =o 


fT 
Moreover, by system equations (4), since aT BYe = 0 and y)” Ae = 0, 
hence we can minimize the error constants corresponding to trees [7] ) and 
[To],, that are given by 


Blo — 2 ZMef’ = (4 - (aT Be) + (oT Be)’ + (aT BMe)”) h8 
= 3 = (a? Be) + (a? Be)’) he, 
2 2 
E{Ioi — 2)" Ze} =(3- ("Ae - (1 Ae) + (Ae) ) h3 


= (4- (7" Ae) + (17 Ae): 3. 


These error constants are minimized with the minimum value 5 if 


1 


at BAe = 5 y?” Ae = 


or equivalently, if 


(2) (3) (3) (6) 


2 
Yi G11 + V2 “Aa + Yq G22 = 9- 


i 


orb!) + azbyy) + aby) = 5, 


By augmenting equations (6) to system (5) and solving the new system by 
MAPLE it is observed that the new system has a three parameters solution 
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that are given by 


a1 = 5 @21 = 59 — 422, A= = 5 
1 1 ao?) 1 1 2) 1 1 

be = a oy Sas ie =- mH, mn ) = —f ), (7) 
2) _ 1” 4@) (8) _@) 2a a 

bit =) b31 = 5 — by: Ja, “= “Ya =3 V2 #0. 


In order to determine the free parameter of the deterministic part, i.e. a22, 
we choose the deterministic part of SRK method (2) to be the Runge —Kutta 
method given by 


that is, it has order 2 (see [4]). This ensures that the semi-implicit method 
works well in the case of small stochastic influence. From (7) we can assume 
A = B®), and consequently for mf) # 0, a one-parameter solution can be 
represented by the following tableau 


1h 0 th +4F, 0 
) 2a tn | Vat in WP Va+in 


1 oe 
In order to choose af ), one can use the minimum of the error constants 


corresponding to trees [[7;],], and [71,71],, that are given by 


T pqy2_72 34.967) 4560461)" 
Elia — 20° Ze] = Sg as h3, 
30727§ 


- oy) Tv 2 (1)2 (1)4 
Elf + $101 — $2 (ZMe)" — $2" Ze] = (=e ) he. 
2 


By introducing two functions f and g in the following form 


(y = 3962? + 560.4 | ) = 148)? +480! 
7 3072 I) = 19944 


it can be shown that these are decreasing functions on the interval (0, +00), 
and moreover 


3 
li a li ae. 
slim, FO) = F95° er Ga 


Now by choosing a) = 3, this class of methods can be represented by the 


following tableau 
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ti. te) 0 
th in| —vB+in —vE+14, 
| th th) -38Vh+h 8Vh 4h 


which is named ‘SJM]2’, and has principal error constants 


tg Voy. Uba0Gs, 39015 
12)’ 12)? 82944’ 82944 


Note that the principal error coefficients corresponding to trees [[7;],], and 
[71,71],, are very close to the limits of f and g as X + +00, respectively. 
Since f and g are even functions, the above analysis shows that the choice 
nf) = 3 is suitable. If we use the 1-norm to estimate the contribution of all 
error terms to the principal error term, then, Tablel represents the following 


values for methods ‘TRK’ ,‘EM1’ , ‘EM2’ and ‘SITM2’. 


Table 1: 1-norm of principal error coefficients 


IRK EM1 EM2 SIM2 
||principal error||; | 1.2083 | 0.89583 | 0.89583 | 0.40019 


From Table1, it follows that the 1-norm principal error of the method ‘SIM2’ 
is less than the 1-norm principal error of ‘EM’ and ‘EM2’ methods. In order 
to improve the results of employing the ‘S/M2’ method at each step, we can 
solve the system for stage-variables Y; and Y2 by the fixed-point iteration 
scheme with starting values for these variables coming from the ‘EM1’ or 
‘EM2’ methods. In fact, for the stage-variable Y; in the ‘SIM2’ method let 


GV) Sm + 5h 90K) + ay (VR+ 6h) (1), 
and hence the fixed-point iteration for solving Y; is given by 
y,8+4 — ayy"), 5s =0,1,2,..., (8) 
with stopping criteria 
yet — Yi" <e, (9) 


where € is a positive known tolerance value. In order to consider the conver- 
gence property of fixed point iterations (8), it is sufficient to have 


GLY) II 5h ob(¥) + SVR + 6h) g(¥) [<1 


Also for the stage-variable Y, let 
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G2(Yo) = Yn + Gh (g0(¥i**4) + go(¥e)) + sie 6A) (a (Mi) + 91 (¥)), 


such that yer satisfy condition (9). Consequently the fixed-point iteration 
for solving Y2 is given by 


¥%"4 = Go(¥%4), 4=0,1,2,.:., (10) 
with stopping criteria 
eer ae ee (11) 


Note that iterations (10) is convergent if 


GAY) [=] Gh ob(¥) + 5o(-Vi + GAY) [<1 


Finally y,+1 for the ‘S1M2’ method will be evaluated by 
1 [s+1] [t+1] 1 [s+1] 
Yne1 =n + 5h (Go (v ) + go (3 )) + —3Vh+ a1 gq (¥! ) 
1 
+ (svi+ =) J (yt) , 


where ylers and yer satisfy conditions (9) and (11). 


3 Numerical results and conclusion 


In this section, the numerical results from the implementation of the above 
seven methods are compared. These methods are ‘IRK’, ‘Milstein’, ‘EM1’, 
‘EM2’,‘SIM1’, ‘IM’ and ‘SITM2’. They will be implemented with constant 
stepsize on two problems taken from [5], for which the exact solution in terms 
of a Wiener process is known. Since J; ~ N(0,h), hence for generating the 
Wiener increments J; in MATLAB environment of random numbers genera- 
tor randn (#traj, #:step) is used, such that each call to randn (#traj, #step) 
creates a #traj x #¢step matrix of independent N(0,1) samples. When these 
methods are simulated, the same sequence of random numbers for the Wiener 
increment J; are used for the stepsize under consideration. The average error 
for each stepsize at the end of the interval of integration is defined by 


K 
1 i : 
AE ==> ly —y (tw) |; 
i=l 
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where yo) is the numerical approximation and y“ (ty) is the exact solution of 
SDE at ty in the 7-th simulation over all K simulations. All of the numerical 
results are based on 1000 simulated trajectories. The results appear in Tables 
2-4. 

Test problem 1. Consider 


dy(t) = —a7y(t)(1 — y*(t))dt + a(1—y*(t))dW(t), y(0) =0, te [0,1] 
with the exact solution 
y(t) = tan h(aW(t) + arctan h(yo)). 
This problem is solved numerically with the choice of parameter a = 1. 


Table 2: Global errors for Test problem 1, with a= 1, K = 1000 and ¢ = 0.001 


T T T T T 
h 35 50 Too 200 700 
TRK 0.21400e-1 | 0.10299e-1 | 0.51948e-2 | 0.24299e-2 | 0.12254e-2 


Milstein |) 0.16276e-1 | 0.82454e-2 | 0.42156e-2 | 0.19930e-2 | 0.10127e-2 
EM1 0.12121e-1 | 0.59344e-2 | 0.30475e-2 | 0.14587e-2 | 0.70585e-3 
EM2 0.12043e-1 | 0.57056e-2 | 0.29270e-2 | 0.13901e-2 | 0.71060¢e-3 
SIM1 0.55857e-2 | 0.21190e-2 | 0.96207e-3 | 0.45136e-3 | 0.22157e-3 

IM 0.13035e-3 | 0.64121e-4 | 0.34962e-4 | 0.17710e-4 | 0.81462e-5 
SIM2 0.80715e-4 | 0.44013e-4 | 0.21736e-4 | 0.10551e-4 | 0.51995e-5 


Test problem 2. Consider 
dy(t) = —(a+ 8 y(t))(1—y* (t) dt + B(1—y?(t))dW(t), y(0) =0.5, te [0,1] 


with the exact solution 


(1) — C+ wo) exp(—2at + 26) + 40 = 1 
WN ~ (1+ yo) exp(—2at + 2BW(t)) +1 — yo 


This problem is solved numerically with a = —1 and for 6 = 1 and 0.01. 
Comparing the numerical results in Tables 2-4, it follows that the ‘S7M2’ 
method is more accurate than the ‘EM1’,‘EM2’, ‘SIM and ‘IM’ methods. 
Also for problems in which the deterministic term dominates (Test problem 
2 with 3 = 0.01) the improvement of the ‘S7M2’ method becomes noticeable 
as the stepsize is reduced. This is because the deterministic component of 
the ‘STM2’ method is the second order Runge-Kutta method. On the other 
hand, for problems in which deterministic term dominates (Test problem 2 
with 8 = 0.01) the global errors for two-stage explicit methods are the same. 
This is because these methods the deterministic components are the same. 
The future work should be based on the construction of implicit IRKs for 
SDEs with two or more Wiener processes. 
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Table 3: Global errors for test problem 2, with a = —1, 8 =1, K = 1000 and e = 0.001. 
h 2 50 rw 20 70 

IRK 0.12763e-1 | 0.58682e-2 | 0.29961e-2 | 0.15034e-2 | 0.74495e-3 
Milstein || 0.11518e-1 | 0.51633e-2 | 0.27770e-2 | 0.13806e-2 | 0.68995e-3 
EM1 0.96413e-2 | 0.41781le-2 | 0.21225e-2 | 0.10660e-2 | 0.54324e-3 

EM2 0.93988e-2 | 0.42298e-2 | 0.20985e-2 | 0.10210e-2 | 0.52317e-3 
SIM1 0.65238-3 | 0.32108e-3 | 0.15186-3 | 0.65537e-4 | 0.30367e-4 

IM 0.79517e-4 | 0.42130e-4 | 0.21167e-4 | 0.10561e-4 | 0.51995e-5 

SIM2 0.57845e-4 | 0.30499e-4 | 0.15806e-4 | 0.79504e-5 | 0.387761e-5 


Table 4: Globa’ 


errors for test problem 2, with a = —1, 8 = 0.01, K = 1000 and e = 0.001 


h 35 a0 ra 200 00 

TRK 0.50778¢e-2 | 0.25193e-2 | 0.12544e-2 | 0.62592e-3 | 0.31264e-3 
Milstein || 0.50778¢e-2 | 0.25193e-2 | 0.12544e-2 | 0.62592e-3 | 0.31264e-3 
EM1 0.50778¢e-2 | 0.25193e-2 | 0.12544e-2 | 0.62592e-3 | 0.31264e-3 

EM2 0.50778e-2 | 0.25193e-2 | 0.12544e-2 | 0.62592e-3 | 0.31264e-3 
SIM1 0.70238e-5 | 0.34182e-5 | 0.17423e-5 | 0.52895e-6 | 0.25624e-6 

IM 0.55102e-5 | 0.26426e-5 | 0.13193e-5 | 0.58242e-6 | 0.29121e-6 

SIM2 0.62103e-6 | 0.15741e-6 | 0.40441e-7 | 0.10595e-7 | 0.28838e-8 
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